Excitonic and Quasiparticle Life Time Effects on Silicon Electron Energy Loss 

Spectrum from First Principles 
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The quasiparticle decays due to electron-electron interaction in silicon are studied by means 
of first-principles all-electron GW approximation. The spectral function as well as the dominant 
relaxation mechanisms giving rise to the finite life time of quasiparticles are analyzed. It is then 
shown that these life times and quasiparticle energies can be used to compute the complex dielectric 
function including many-body effects without resorting to empirical broadening to mimic the decay 
of excited states. This method is applied for the computation of the electron energy loss spectrum 
of silicon. The location and line shape of the plasmon peak are discussed in detail. 
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One of the most challenging condensed matter prob- 
lems is the prediction of excited states of materials from 
first principles lj . The solution of this formidable task is 
doubly rewarding since (1) most of the interesting physics 
involves the interaction of an electromagnetic radiation 
with matter, and (2) most of progress in today's nano- 
electronic technologies requires the knowledge of excited 
states. 

In the last few decades, the development of this field 
has been expending rapidly, and several evolutionary pe- 
riods can be distinguished. In the first one, the den- 
sity functional theory in the local density approximation 
(LDA) played a major role, and was widely used for the 
analysis of spectroscopic properties of materials 0,0, de- 
spite that the band gaps of semiconductors and insula- 
tors were adjusted by ad hoc methods. Later, it was 
realized that the GW approximation of HedinQ pro- 
vides a practical scheme for the determination of the 
quasiparticle energies, and band gaps were no longer 
adjusted 5]. Nevertheless, even the use of correct band 
gaps and the various types of exchange-correlation cor- 
rections, including local-field effects, did not improve sig- 
nificantly the calculated optical spectraQ. It was often 
believed that the inclusion of the electron- hole ( e-h) in- 
teractions was the missing ingredient for an adequate de- 
scription of the optical spectra, and this was later con- 
firmed by model calculations 0. However, it is only re- 
cently that ab initio methods explicitly showed the rel- 
evance of these complex interactions to the calculation 
of the dielectric function pj. At this point, calculations 
from first-principles were finally making direct contact 
with experiment and establishing the relevance of the e- 
h interactions. 

The aim of this Letter is to add a further step to 
this fascinating development by incorporating, from first 
principles and for the first time, the life time of inter- 



band transitions into the calculation of the dielectric 
function. As a consequence an empirical broadening is 
no longer needed for a successful comparison of theoret- 
ical and experimental optical spectra. This work also 
improves upon existing excited state calculations based 
on the pseudopotential (PP) approach in conjunction 
with plasmon-pole (PIP) models for the screening of the 
Coulomb interaction^] . In those types of calculations 
not only the various types of matrix elements of operators 
are not accurate due to the use of pseudowave functions, 
the PIP approximation also makes it impossible to deter- 
mine the imaginary part of the self-energy, and hence the 
spectral functions and life times of quasiparticles remain 
unaccessible. Indeed, it has been noticed recentlyHHI3 
that GWA implementations based on PP methods led to 
larger and more k-dependent shifts than those based on 
all-electron methods. To illustrate our method, we com- 
pute the plasmon resonance and the line shape of the 
electron energy loss spectrum (EELS) of silicon, one of 
the most used materials in today's technology. 

The most successful approach for describing elemen- 
tary excitations probed by photoemission experiments 
uses many-body Green's function theory. The energies 
and life times of quasiparticles are mainly determined by 
the pole structure of the Green's function G or more con- 
veniently by solving the quasiparticle equation^}]. The 
self-energy operator E = EW + iE^ 2 ' of the quasiparti- 
cles is non-local, frequency-dependant, and generally non 
hermitian. The non-hermitian part iE^ 2 ) is related to the 
dominant relaxation mechanisms (scattering on crystal 
imperfection, electron-phonon or electron-electron inter- 
action) which give rise to the finite life time of quasi- 
particles. In the present work, E is computed within 
Hedin's GW framework 4] which includes dynamic po- 
larization in the random-phase approximation (RPA). In 
such an approach, only the electron-electron interaction 
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FIG. 1: Spectral function of silicon for L, V and X k points. 
The q vectors are given in units of 27r/a, where a =10.261 
a.u. The energy zero is fixed at the energy of the topmost 
occupied LDA state at the V point. 



contributes to the non-hermitian part of the self-energy, 
i.c to the instability of single-particle excitation. 

The complex QP energies e?£ are solutions of the 
quasiparticle equation for a state labeled nk, and are 
in practice obtained using a first order perturbation 
theory 0. The real part of corresponds to the QP 
energy, and the life time of a single-particle excitation, 
given by the inverse of the full-width at half-maximum 
(FWHM) of the QP peak in the spectral function, is 
defined by r„ k = [2 x |r„ k |] _1 where r„ k = Im[e^]. 
The determination of both QP energies and their life 
times requires the computation of the self-energy matrix 
elements [ri|. 

The LDA eigenvalues and eigenvectors obtained by 
means of the all-electron Projector Augmented- Wave 
method^3 are used as a starting point for our GWA cal- 
culations. The wave functions are expanded into plane- 
waves up to 20 Ry and three partial waves of s, p and d 
types are used to describe the correct nodal structure of 
the wave functions near the nuclei. 

In Fig. ^the spectral function of some selected high 
symmetry k-points (L, T, X) of silicon is displayed. The 
FWHM of the QP peaks increases almost linearly away 
from the band gap because more and more decay chan- 
nels of Auger type become available. Physically speaking, 
a QP is scattered into LDA empty states whose energies 
are between the QP energy and the Fermi energy. 

The GW self-consistency should affect the QP life time. 
Indeed, the number of decay channels of Auger type is 
fixed by the smallest possible excitation energy contained 
in the screened interaction W(u>) which is defined by the 
absolute band gap of silicon. As the band gap increases 
from 0.44 eV to 1.0 eV when going from the LDA to 
the GW approximation (see Fig. 0a), we expect a de- 
crease of the imaginary part of the self-energy caused by 
self-consistent effects. Therefore the QP peak should be 
narrowed, corresponding to increased life times. Another 
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FIG. 2: (a) Shifts of QP peak positions Re[e qp ] - e LDA of 
silicon as function of the LDA energies e LDA for the first 16 
bands and for 216 k-points sampling the full Brillouin zone, 
(b) Imaginary part |7m(e 9P )| and life times r qv of QPs as a 
function of the QP peak positions. The energy zero is fixed at 
the energy of the topmost occupied LDA state at the P point. 

consequence of the lack of self-consistency is that some 
states near the Fermi level have a spurious imaginary 
part because they fall between the LDA and GW Auger 
thresholds. Thus, the life times of the QPs in the vicin- 
ity of the gap are not correct except at the uppermost 
valence state and at the lowest conduction states where 
the non-hermitian part of the self-energy due to phonons 
is strictly zero 13] . A proper description of the QP life 
time in this range of energy should include the interac- 
tion of electronic states with lattice vibrations because 
the dominant relaxation process is governed by electron- 
phonon interaction. We expect, however, that our predic- 
tions are quantitatively correct at high energy where the 
dominant relaxation mechanism is electron-electron in- 
teraction. Therefore, the complex QP energies obtained 
by means of the GW approximation represent a good 
starting point for the direct calculation of the complex 
dielectric function at high energy without using any phe- 
nomenological life time of excited states. 

To include the e-h interaction in the dielectric function 
we solve the Bethe-Salpeter equation, which describes 
the propagation of a correlated quasielectron-quasihole 
pair, and can be turned into an effective two-particle 
Schrodinger Hamiltonian whose block matrix form is 
given by 

= ( _% _% ) . (1) 

Detailed expressions of the different block matrix ele- 
ments can be found elsewhere [l4|. The diagonal blocks 
given by R and — R* are respectively the resonant part 
(i>k — » ck transitions) and the anti-resonant part (ck — > 
uk transitions) while the off-diagonal blocks C and — C* 
couple positive and negative frequency transitions. 01- 
evano and Reining have shown that these off-diagonal 
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blocks are crucial for the correct description of the EELS 
of silicon [l5|. Here, we follow a similar approach, but in 
contrast to their work, we include the life time of non- 
interacting e-h pairs by considering complex transition 
energies whose imaginary part is set by Im[e^ — e*£] = 



r c k — r„k- Therefore, even if one neglects the off-diagonal 
blocks, the effective Hamiltonian H cxc is non-hermitian 
because the diagonal matrix elements of R are complex. 

The macroscopic dielectric function is given by 



e(q = 0, W ) = 1- hm - x x £ (mle^-'K) x [H™ * (n 4 |e iq - r |n 3 ) x (/„ 4 -/„ 3 ), (2) 



with f nt being the occupation number of the Bloch state 
labeled by ni and with (nj,ny) = (vis., ck) or (ck, t>k). 

To build the excitonic Hamiltonian H exc , we first cal- 
culate the complex QP energies for the first 16 bands and 
for the 216 k points sampling the Brillouin zone. Fig. |2Ja 
shows the shifts of the real part of the QP energies with 
respect to the LDA energies as a function of LDA ener- 
gies. Interestingly, these shifts are far from being uniform 
across the BZ with the exception of LDA states ranging 
from -5 eV to 5 eV. In this range of energy, the valence 
LDA states are shifted downward by about 0.45 eV while 
the LDA conduction states are shifted upward by about 
0.15 eV. The striking feature is the decrease of the QP 
energy shifts for the LDA conduction states above 5 eV. 
This trend was also observed by Fleszar and Hanke[l6j. 
In addition, the QP shifts below -5 eV reflect the con- 
traction of the occupied GW bandwidth with respect to 
that of the LDA. Such a contraction is not observed when 
using a plasmon-pole model. Fig. |5Jb shows the imagi- 
nary part of the QP as a function of the QP peak posi- 
tions. It's worth mentioning that the imaginary part of 
the QP of Ref. are twice as larger as ours. Such a 
large difference can certainly be ascribed to the value of 
the broadening parameter S used in Ref. 0] to evaluate 
W(uj) along the real axis (see Ref. for details). The 
broadening should be chosen as small as possible (6 — > 0) 
to obtain correct life times. Notice that |Im(e^)| is very 
small for the QP states ranging from -5 eV to 5 eV and 
increases almost linearly outside this interval. Thus, we 
expect a non negligible broadening |r c k — T^kl of about 
0.5 eV for non interacting transition energies in the range 
of energy where the plasmon resonance occurs. 

Fig. El shows our calculated EELS (— Im[e _1 (q = 
0, uj )]) of silicon[l3. These results are compared with 
the experimental™ectrum[ll characterized by a rather 
symmetric plasmon peak located at 16.7 eV. The dashed 
line represents the spectrum calculated within the RPA 
approximation which amounts to retaining just the com- 
plex transition energies in the excitonic Hamiltonian. 
The main peak of the RPA spectrum is shifted by about 
1 eV towards higher energy compared to experiment. In 
addition, the FWHM of the plasmon peak AEi is un- 
derestimated by a factor of 3 while the height is overesti- 
mated by the same factor. When including both e-h at- 
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FIG. 3: EELS of silicon: the thick line is the experimental 
data by Stiebling (Ref. f!s||1. the dashed line is the RPA cal- 
culation without excitonic effects (local-field effects are also 
neglected), the curve with empty circles includes both exci- 
tonic effects and local-field effects, while the full circle curve is 
same as empty circle one but includes the coupling. The com- 
plex GW eigenvalues have been used in all calculations. The 
real part ei and imaginary part £2 of the dielectric function 
are shown in the inset (Arrows show the zero of ei(w)). 



traction and local-field effects but still neglecting the cou- 
pling between positive and negative frequency transitions 
(open circles), the line shape of the spectrum worsens. 
The plasmon-peak position is now ill-defined because the 
spectrum exhibits two structures located at 15 and 18 
eV. Finally, the inclusion of the coupling between for- 
ward and backward going quasielectron-quasihole pairs 
by taking into account the off-diagonal blocks of H exc 
has a rather drastic effect on the overall spectrum (full 
circles). On the one hand, the natural line shape is im- 
proved and, as expected, AEi ~ 3 eV is underestimated 
with respect to experiment (~ 3.7 eV) because of the 
analyzer resolution included in the latter. On the other 
hand, the main peak of our spectrum is unfortunately 
shifted, by about 1 eV, towards lower energy compared 
to experiment, and its height is somehow overestimated. 
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However, the inclusion of the analyzer resolution would 
slightly reduce the calculated height. 

The position of the plasmon peak lu p is determined 
roughly by the zero of the real part e\ (w) of the macro- 
scopic dielectric function. The inset in Fig. [3] clearly 
shows the shift of the zero towards lower energies when 
the excitonic effects and coupling terms are included. 
It's also of interest to notice that the FWHM of the 
plasmon peak is qualitatively given by AEi ~ 2 x 

e 2(k>p)/ 8< g^ Up- Thus AEi is controlled by the value 
of £2 at the plasmon resonance energy and by the slope 
of ei at the same energy as it can be seen in Fig. [3] 
For example, a detailed comparison of both calculations, 
where the exchange-correlation effects are taken into ac- 
count, shows that the inclusion of the coupling yields a 
decrease in ti as well as an increase in the slope of e±, 
i.e., it leads to a narrowing of the plasmon peak. 

Obviously, the puzzling feature of our calculated spec- 
trum is the underestimation of the plasmon-peak energy. 
It is therefore worth comparing our findings to the only 
available PP results[L5|. The PP plasmon peak is located 
at about 17.3 eV while ours is located at 15.7 eV. These 
differences arise because we avoid most of the approx- 
imations used in Ref. fl5l In particular, (1) we avoid 
using PP methods, and any kind of plasmon-pole model 
to calculate QP energies (such models are expected to fail 
at these high energies), (2) we don't use any broadening 
parameter to describe the life times of excited states, and 



finally (3) we avoid using first order perturbation theory 
to include the coupling between positive and negative fre- 
quency transitions. This latter approximation is shown 
to have only a minor effect on the final spectrum |15| . 
As a consequence of the first two approximations, their 
plasmon-peak energy position, at the RPA level with- 
out excitonic effects, is already shifted by 0.9 eV towards 
higher energies with respect to ours. A self-consistent 
GW calculation should improve the agreement with 
experiment by enlarging the transition energies. How- 
ever, self-consistency issues are beyond the scope of this 
work. 

In conclusion, we have calculated the QP energies and 
life times of silicon within the all-electron GW approxi- 
mation. We have also demonstrated that these QP ener- 
gies can be used to evaluate the EELS of silicon without 
using any phenomenological parameter to describe the 
life times of excited states. Our results confirm that it 
is crucial to include excitonic effects and coupling be- 
tween forward and backward going e-h pairs to obtain 
an improved line shape. However, the underestimation 
of the plasmon-peak position reflects the need for a self- 
consistent calculation of the QP energies or the need for 
an approximation beyond the GW method. 
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time has been granted by CINES (Project No. gmg2309) 
on the IBM SP4. 
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